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Counterions at Charged Cylinders: Criticality and universality beyond mean-field 
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The counterion-condensation transition at charged cylinders is studied using Monte-Carlo sim- 
ulation methods. Employing logarithmically rescaled radial coordinates, large system sizes are 
tractable and the critical behavior is determined by a combined finite-size and finite-ion-number 
analysis. Critical counterion localization exponents are introduced and found to be in accord with 
mean-field theory both in 2 and 3 dimensions. In 3D the heat capacity shows a universal jump at the 
transition, while in 2D, it consists of discrete peaks where single counterions successively condense. 
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Many biopolymers, such as DNA, actin, tubulin, fd- 
viruses, are charged and stiff. On length scales smaller 
than the persistence length, they can be represented by 
straight, charged cylinders and oppositely charged ions 
(counterions) are attracted via an electrostatic potential 
that grows logarithmically with radial distance. As the 
ion confinement entropy also exhibits a logarithmic de- 
pendence, it was suggested early by Onsager that a coun- 
terion derealization transition occurs at a critical cylin- 
der charge (or equivalent at a critical temperature) Q . 
This argument is strictly valid only for a single ion since it 
neglects cooperativity due to inter-ionic repulsions. Nev- 
ertheless, it was corroborated by mean-field (MF) studies 
HI El , which demonstrate that below a critical tem- 
perature, a fraction of counterions stays bound or con- 
densed in the vicinity of the central cylinder even in the 
limit of infinite system size; while above the critical tem- 
perature, all counterions de-condense to infinity. This 
counterion-condensation transition (CCT) dramatically 
affects a whole number of static and dynamic quantities 
for charged polymers . It has been observed with differ- 
ent polymers by varying the medium dielectric constant 
0, |i| or the polymer charge density 0, |(| ; the counte- 
rion distribution around DNA strands has been directly 
measured recently using anomalous scattering techniques 
• Since its discovery, the CCT has been at the focus of 
numerical 8, j| and analytical ^} studies. Under partic- 
ular dispute has been the connection between CCT and 
the celebrated Kosterlitz-Thouless transition of logarith- 
mically interacting particles in 2 dimensions Also, 
the precise location of the CCT critical point remains 
subject of ongoing experimental investigations 

As is well known from bulk critical phenomena, fluc- 
tuations and correlations typically make non-universal 
and universal quantities deviate from mean-field theory 
(MFT) below the upper critical dimension [12|. Surpris- 
ingly, the MFT prediction for the CCT critical tempera- 
ture has not been questioned in literature and apparently 
assumed to be exact. Likewise, the existence of scaling 
relations and critical exponents associated with the CCT 
has not been considered, neither on the MFT level and 
(consequently) also not in the presence of correlations. 



In this paper we pose the questions: i) what is the 
critical temperature of the CCT, and ii) what are the as- 
sociated relevant critical exponents? We employ Monte- 
Carlo simulations, which are performed in rescaled log- 
arithmic coordinates in order to handle very large sys- 
tems (where the criticality actually occurs) with tractable 
equilibration times. A combined finite-size and finite-ion- 
number analysis yields the desired critical temperature 
and exponents. To enhance the effects of fluctuations, we 
also study a (within MFT equivalent) 2D system of loga- 
rithmically interacting charges, as applicable to an exper- 
imental system of oriented cationic and anionic polymers 
(e.g. DNA with polylysine 13] ). Surprisingly, MFT is 
demonstrated to be accurate both in 3D (charged cylin- 
der with point-like counterions) and in 2D (charged cylin- 
der with cylindrical counterions). The critical exponents 
associated with the inverse counterion localization length 
(which plays the order parameter of the CCT) and the 
universal behavior of the heat capacity are determined; 
both quantities are experimentally accessible. 

In the 3D simulations, we consider a central cylin- 
der of radius R and uniform surface charge density a s 
(linear charge density r = 2irRa s ) with N neutraliz- 
ing point-like counterions of valency q confined laterally 
in an outer cylindrical box of radius D (see Fig. ^ for 
snapshots projected along the z-axis). Periodic bound- 
ary conditions in ^-direction are handled using summa- 
tion methods for long-range interactions [l4|. Rescaling 
all spatial coordinates by the Gouy-Chapman length, 
fj, = l/(27r<7^BC r s), as x = x/ fi, we obtain the Hamil- 
tonian H = 2^J2f=i Mn/R) +^J2{ij) Vl** ~ %l (P er 
k-QT), where £b = e 2 / (ATreeoksT) is the Bjerrum length, 
and f is the radial distance from the cylinder axis. 
The coupling parameter, 5 = 2irq 3 £ 2 i o- s , is an indica- 
tor for the importance of ionic correlations: in the limit 
S — > 0, correlations are unimportant and MFT becomes 
exact; in the converse strong-coupling limit 5 — * oo, 
MFT breaks down [l^. The so-called Manning pa- 
rameter (rescaled inverse temperature), £ = ql^r, reg- 
ulates the CCT and is a measure of counterion bind- 
ing: According to MFT 0, 0, 0, the CCT occurs 
for A = ]n(D/R) — > oo at the MF critical threshold 
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FIG. 1: Snapshot-topviews for a) logarithmic box radius A = \n(D/R) = 100 and Manning parameters £ = 0.7, 1.0 and 3.0, and 
b) £ = 1.0 and A = 10 and 25 shown in logarithmic radial units y = ln(r/7?) (number of particles TV = 100 and the coupling 
parameter S = 0.1). c) Radial counterion distribution function, p(f), for £ = 1.0, S = 0.1, TV = 100 and various system sizes, 
d) Counterion distribution for £ = 3.0 and various 3 shows a crossover between MFT (solid curve) and SC (dashed curve) 
predictions, Eqs. J2J and ©. e) Order parameter Si = (1/r) as a function of Manning parameter, £, for different coupling 3 
compared with MFT (solid curve) and SC (dashed curve) results, Eqs. J3J and JjJ. In d), e) A = 300 and typically TV = 200. 



^mf __ abovg which a fraction 1 — l/£ of counteri- 
ons condenses. To investigate the critical limit for large 
lateral system size, we introduce a (centrifugal) sam- 
pling method by mapping the radial coordinate to the 
logarithmic scale as y — ln(r/R). The partition func- 
tion transforms as Z ~ J fjdfjdfjd^i] exp(— TL) ~ 
/ Uli d?/idzid0i] exp(— Kmc), where the "Hamiltonian" 



Hmc — 2(£ 



AT 

i=l 



(1) 



is used for Monte-Carlo sampling; it features a linear po- 
tential (first term) acting on counterions from competing 
energetic (~ 2£y) and entropic or centrifugal (~ 2y) con- 
tributions associated with the cylindrical boundary. 

In Fig. ^i, we show the snapshot-topviews from our 
simulations for A = \n(D/R) = 100. De-condensation 
phase is reproduced for small Manning parameter, £ — 
0.7, as counterions gather at the outer boundary; while 
for £ = 3, a fraction of counterions accumulates or con- 
denses around the central cylinder. The transition regime 
for intermediate £ exhibits strong finite-size effects: As 
seen for £ = 1 in Figs.^, and b, only for large logarithmic 
system size, A ^> 1, does de-condensation occur. This is 
also demonstrated by vanishing radial distribution func- 
tion of counterions, p(f), for growing A in Fig. ^j. For 
small coupling parameter 2 = 0.1, the data for p[r) com- 
pare well with the normalized MFT profile (solid curves) 
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where £ > A/(l + A) and (3 is given by £ = (1 + /3 2 )/(l - 
/3cot(— /3A)) 0| . Conversely, for large coupling 3 — > oo, 



strong-coupling (SC) theory becomes valid and yields 
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The break-down of MFT for increasing coupling param- 
eter is demonstrated in Fig. (at fixed £ = 3), where 
the distribution function exhibits a gradual crossover be- 
tween the asymptotic MFT and SC predictions [l5| . En- 
hanced ionic correlations at elevated 2 cause a remark- 
ably larger counterion density near the charged cylinder 
■ T ne central question is whether these correlations 
influence the critical behavior associated with the CCT. 

The standard analysis of critical behavior relies on 
suitably-defined order parameters 0]. Here, this role 
is taken by the inverse moments of the counterion distri- 
bution function, S n — (1/r") (with n > 0); S n vanishes 
in the de-condensation phase, but attains a finite value in 
the condensation phase. The MFT distribution J2J) gives 
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while in the strong-coupling limit, we have from Eq. J3J 

sc = 2(g-l) l- e -(2€-2+n)A 

e(2e-2 + n) l-e- 2 («- 1 ) A ' 



(4) 

3 

(5) 



The data for S\ (the mean inverse localization length) 
in Fig. ^ exhibit the condensation (Si > 0) and de- 
condensation (Si = 0) phases for a wide range of cou- 
plings, 2. To locate the critical Manning parameter, £ c , 
a finite-size analysis is required since criticality is masked 
both by finite system size, A, and finite particle number, 
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FIG. 2: a) Rescaled (potential) energy, E / A, and b) rescaled (excess) heat capacity, C/A, as a function of Manning parameter, 
£, for fixed N = 100 compared with the MFT prediction for A — > oo (solid line), A = 300 and 100 (broken lines, from top to 
bottom) [T^| . c) Location of the peak of the energy, £jf, as a function of 1/A for S = 0.1 compared with MFT (solid line); 
Inset: Peak location versus S for TV = 200 and A = 300. d) Scale invariance of the order parameter, Si, near £ c with respect 
to the reduced Manning parameter, f = 1 — £ c /£, and the counterion number, N, with exponents a/c = 2/3 and 1/c = 1/3 
(lateral system size A = 300). e) Scale invariance of Si for various A with exponents a/b = 2 and 1/6 = 1 (N = 200 is fixed). 



N, in the simulations (i.e. S\ does not vanish and sat- 
urates at a small value at critical point). To this end, 
we study the singular behavior of rescaled energy E — 
{H/(Nk B T)) and heat capacity C = {(SH / k B T) 2 ) / N , 
where 8TL = H — (TL). Simulation results for E and C 
(Figs. |2K, b) show a non-monotonic behavior that can 
be understood for large systems, A » 1, using a simple 
argument: for small £, all counterions are unbound and 
the electrostatic potential experienced by counterions is 
that of a bare cylinder, ip(f) — 2£ln(f/i?) (per k B T/eq), 
at the outer boundary. The energy follows by a charging 
process as E ~ ip(D)/2 ~ £A. From the thermody- 
namic relation £dE/dt; = E — C, one obtains C ~ 0. For 
large £, the cylinder potential is screened by condensed 
counterions. Using the MFT relation p ~ e"^, Eq. (J2J 
gives ip(f) ~ 2 Inf. A fraction l/£ of counterions is de- 
condensed 0,0, E| an d yields the leading contribution to 
the energy; hence E ~ ip(D)/(2g) ~ A/£ and C ~ 2A/£. 
Thus both quantities decay with £ as counterions become 
increasingly condensed. The energy exhibits a peak and 
the heat capacity develops a jump, which become singu- 
lar for an infinite system (A — > oo, N — > oo) reflecting 
the CCT point, £ c (Figs. b). We determine £ c from 
the location of the energy peak, £f (A, N), for increasing 
N and A. The numerical results for £f (A, N) in Fig.|2t, 
obtained using the identity £<9-E/9£ = E — C, compare 
favorably with the MF prediction (solid curve) for small 
coupling 5 = 0.1 and increasing N. The MF prediction 
asymptotically tends to the MF threshold £^* F — 1 ac- 



jump) shows no dependence on the coupling parameter 
within error bars (Fig.[5J; inset). We thus find a universal 
counterion-condensation threshold as £ c = 1.00 ± 0.002. 

We now turn to the near-threshold scaling behavior 
of the order parameter. In the thermodynamic infinite- 
dilution limit (N — ► oo, A — > oo), S n exhibits a power- 
law behavior as S n ~ ( x (close to and above £ c ), where 
£ = 1 — £ c /£ is the reduced Manning parameter (reduced 
temperature). Within MFT, we obtain xmf = 2 for all 
n from Eq. Q. At criticality, £ — 0, one expects the 
scaling S n ~ A" 7 for N -> oo but finite A (within MFT, 
we obtain 7mf = 2): while for A — > oo but finite N, one 
expects S n ~ N~ v (v is not defined in MFT). These rela- 
tions indicate that close to criticality, S n ((, A, N) takes 
a homogeneous scale-invariant form, i.e. for A > 0, 



S n (A£, A~ b A, \- c N) = \ a S n ((, A, N) 



(6) 



cording to £ 



MF 



£f' MF (A) ~ 1/A as A 
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location of the energy peak (and also the heat capacity 



where a, b, c are related to the exponents 7, v, x- Choos- 
ing the scale factor as A = N x / C , one finds S n (£, A, N) = 
7V^ Q / c S' n (CiV 1/c , AN- fc / c , l). For large A7V- fc / c , as is the 
case in our simulations, and at the transition, £ = 0, we 
have S n ~ N~ a / C and thus obtain v = a/c. For N — > 00, 
a similar argument leads to jSV» ~ £ a , which gives x = a - 
In Fig. Eli we plot the rescaled order parameter N a /°S\ 
as a function of (N 1 ^ for various N and fixed large A. 
By choosing the scaling exponents as v = a/c — 2/3±0.1 
and x — a = 2 ± 0.4, we obtain excellent data collapse 
both for small coupling 5 (main figure) and large 5 (in- 
set). Choosing A = A 1 ^ in Eq. © yields the relation 
S n ((,A,N) = A- a / b S n ((A 1 / b ,l,NA- c / b ). InFig.^we 
plot A a / b Si versus C t A 1 / h . Good data collapse is obtained 
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FIG. 3: 2D results: a) Order parameter Si = (1/r) as a func- 
tion of Manning parameter, £, for different particle numbers, 
TV, compared with MFT and SC predictions, Eqs. Q) and JSJ. 
b) Rescaled energy and c) heat capacity (per particle) exhibit 
a set of singularities corresponding to successive condensation 
of single counterions for increasing f . Here A = 300. 

for 7 = a/b = 2 ± 0.6 both at small (main figure) and 
large 5 (inset), which also demonstrates approximate in- 
dependence from the scaling argument NA~ c / b 16]. Our 
numerical results give the same critical exponents, 7, v 
and x> f° r a U n an d agree with MFT predictions (for 
7 and x)- The exponents appear to be universal, i.e., 
independent of the coupling parameter, 3. 

Deviations from MFT in general grow with diminish- 
ing dimension 12]. The two-dimensional analogue of 
the counterion-cylinder system consists of logarithmically 
interacting mobile counterions and a central charged 
disk. The 2D Hamiltonian reads Hn — 2£X}i=i — 
ln(|xj — Xj|). Unlike in 3D, £ and 5 are related 
due to electroneutrality as 3 = £/N. Thus the strik- 
ing feature in 2D is that for a given Manning parameter, 
£, the coupling parameter tends to zero, 5 — > 0, in the 
limit of many counterions, N — > 00, and MFT should 
become exact. Fig. [3^ shows the simulated order param- 
eter Si in 2D. For N — 1, the data trivially follow the 
SC prediction dashed curve, and for increasing N, 
they tend to the MF prediction (J2J, solid curve. Accord- 
ingly, scaling analysis of the condensation threshold and 
the critical exponents for N — ► 00 gives identical results 
as in 3D and thus no deviations from MFT. However, 
closer inspection of the 2D data in Fig. reveals a pe- 
culiar set of cusp-like singularities for finite N. In fact, 
these singularities correspond to delocalization events of 
individual counterions, which give rise to a sawtooth- 
like structure for mean energy and a series of discrete 



peaks for heat capacity (Figs. Ob, c). This can be under- 
stood by a simple analysis of the 2D partition function: 
Suppose that N — m counterions are firmly bound to 
the central cylinder (disk), while m — 1 counterions have 
evaporated to infinity (m = l,...,iV). Neglecting the 
delocalized ions, the partition function can be written 



as Z N = J Ili 



=m+l 



exp(-7YAr„ m ) x Z N 



< ' m \ where 
l ln(|xj— x])] is the 



contribution from the m-th counterion which is assumed 
to be weakly localized. It is thus de-correlated from 
the firmly bound ions and 2^ approximately factorizes 



In 1x1' 



jy diverges for 
For N counterions this 



00, Z\ 



as Z^P' — J d aiexp[— 2^1ni 
e (2-2me/iV)A > In the Umit A . 

Manning parameters ^ < N/m 
gives a discrete set of singularities at £ m = N / m with di- 
verging heat capacity C ~ (^ — ^m)^ 2 , and in agreement 
with our simulations (Fig. |3J)- For a Manning parameter 
range N/(m + 1) < £ < N/m, there are m de-condensed 
ions. In the thermodynamic limit N — > 00, the fraction 
of de-condensed ions, m/N, becomes a smooth function 
and tends to the MF prediction jHQ, i.e., m/N — > l/£. 

In summary, both in 2D and 3D the location and crit- 
ical exponents of the counterion-condensation transition 
at charged cylinders are correctly described by mean-field 
theory. The heat capacity is experimentally accessible: in 
2D (parallel charged polymers), it consists of a discrete 
set of peaks at which single counterions condense, and in 
the limit N — > 00, it converges to the 3D shape with a 
universal jump at the condensation threshold. 

Discussions with H. Boroudjerdi, Y. Burak, H. Orland, 
Y. Y. Suzuki and support from the DFG (SFB-486) are 
acknowledged. 
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